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A  novel  parametrized  model  that  describes  the  3D  microstructure  of  compressed  fiber-based  materials  is 
introduced.  It  allows  to  virtually  generate  the  microstructure  of  realistically  compressed  gas-diffusion 
layers  (GDL).  Given  the  input  of  a  3D  microstructure  of  some  fiber-based  material,  the  model  com¬ 
presses  the  system  of  fibers  in  a  uniaxial  direction  for  arbitrary  compression  rates.  The  basic  idea  is  to 
translate  the  fibers  in  the  direction  of  compression  according  to  a  vector  field  which  depends  on  the  rate 
of  compression  and  on  the  locations  of  fibers  within  the  material.  In  order  to  apply  the  model  to 
experimental  3D  image  data  of  fiber-based  materials  given  for  several  compression  states,  an  optimal 
vector  field  is  estimated  by  simulated  annealing.  The  model  is  applied  to  3D  image  data  of  non-woven 
GDL  in  PEMFC  gained  by  synchrotron  tomography  for  different  compression  rates.  The  compression 
model  is  validated  by  comparing  structural  characteristics  computed  for  experimentally  compressed  and 
virtually  compressed  microstructures,  where  two  kinds  of  compression  -  using  a  flat  stamp  and  a  stamp 
with  a  flow-field  profile  -  are  applied.  For  both  stamps  types,  a  good  agreement  is  found.  Furthermore, 
the  compression  model  is  combined  with  a  stochastic  3D  microstructure  model  for  uncompressed  fiber- 
based  materials.  This  allows  to  efficiently  generate  compressed  fiber-based  microstructures  in  arbitrary 
volumes. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

We  introduce  a  novel  parametrized  model  that  describes  the  3D 
microstructure  of  compressed  fiber-based  materials.  This  model  is 
of  high  practical  relevance  since  it  allows  for  generating  the 
microstructure  of  gas-diffusion  layers  (GDL)  under  consideration  of 
compression  conditions  as  found  in  fuel  cells.  Given  the  input  of  a 
3D  microstructure  of  some  (uncompressed)  fiber-system,  the 
model  compresses  the  system  of  fibers  in  a  uniaxial  direction.  By 
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the  term  ‘uniaxial  compression’  we  mean  that  the  material,  located 
on  a  fixed  ground,  is  compressed  by  a  stamp  which  presses  the 
material  from  a  uniaxial  direction  against  the  fixed  ground,  see 
Fig.  1.  The  compression  rate  is  then  given  by  1  minus  the  ratio  of  the 
thickness  of  the  compressed  material  and  the  thickness  of  the 
uncompressed  material.  The  model  can  be  applied  to  mimic  any 
uniaxial  compression  of  fiber-based  materials.  In  particular,  it  is 
able  to  describe  the  compression  for  varying  stamp  geometries  and 
for  arbitrary  compression  rates.  In  addition,  when  applied  to  sys¬ 
tems  of  non-overlapping  fibers,  i.e.  fibers  which  do  not  penetrate 
each  other,  the  compression  model  keeps  the  fiber  system  free  of 
overlaps  and  it  preserves  the  total  volume  of  the  fibers.  Note  that 
the  idea  of  the  compression  model  introduced  in  the  present  paper 
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Fig.  1.  3D  image  of  non- woven  GDL  (yellow  fibers)  with  flat  stamp  (left)  and  with  stamp  of  flow-field  structure  (right). 


is  to  rebuild  the  structure  of  compressed  materials,  i.e.,  it  is  a  purely 
morphological  approach  which  does  not  include  any  physics.  The 
advantage  of  this  morphological  approach  is  that  compression  of 
materials  can  be  mimicked  by  means  of  simulations  in  short  time 
on  huge  domains  and  additionally  the  simulations  based  on  the 
compression  model  are  good  approximations  of  real  compressed 
materials. 

We  apply  the  compression  model  to  3D  image  data  describing 
the  microstructure  of  non-woven  GDL  in  proton  exchange  mem¬ 
brane  fuel  cells  (PEMFC),  which  consists  of  (non-overlapping)  fi¬ 
bers.  PEMFC  are  a  seminal  technology  for  electrical  power 
generation  due  to  their  high  efficiency  and  environment  friendly 
emissions  [1,2].  An  important  factor  for  the  efficiency  and  stability 
of  the  PEMFC  is  its  GDL,  which  is  mainly  responsible  for  the  gas 
supply  of  the  electrodes  and  water  storage/evacuation  within  the 
GDL.  For  further  details  concerning  the  functionality  of  GDL  we 
refer  e.g.  to  [3-7].  The  GDL  under  examination  is  a  porous  material 
whose  solid  phase  consists  of  a  system  of  strongly  curved  (non¬ 
overlapping)  carbon  fibers.  The  microstructure  of  the  GDL  is  closely 
related  to  its  functionality,  in  particular  to  the  transportation  of 
gases  and  fluids  through  the  GDL  [8-12].  Thus,  the  systematic 
development  of  ‘designed’  morphologies  with  improved  function¬ 
ality  is  an  important,  but  largely  unsolved  task. 

In  the  operating  state  of  PEMFC,  the  GDL  is  uniaxially  com¬ 
pressed  under  the  lands  of  the  flow-field  in  order  to  ensure  the 
mechanical  support  of  the  membrane  as  well  as  thermal  and 
electrical  conductivity.  It  was  demonstrated  that  the  compression 
rate  of  the  GDL  influences  the  performance  of  a  PEMFC  [13-17  .  In 
Fig.  1  (right),  the  microstructure  of  a  non-woven  GDL  is  displayed 
with  a  stamp  possessing  a  flow-field  structure.  In  order  to  improve 
the  understanding  of  transportation  processes  within  the  GDL  in 
operating  PEMFC,  it  is  of  great  importance  to  investigate  the  cor¬ 
relation  between  compressed  microstructures  of  GDL  and  the  cor¬ 
responding  physical  processes,  e.g.  gas  and  water  transport. 

The  first  step  towards  the  analysis  of  this  correlation  is  to  un¬ 
derstand  how  the  microstructure  changes  under  compression.  In 
particular,  there  is  a  need  of  models  describing  the  3D  micro¬ 
structure  of  compressed  GDL  by  the  flow-field  for  varying  rates  of 
compression. 

In  general,  there  exist  two  classes  of  models  which  virtually 
compress  microstructures.  The  first  class  is  based  on  numerical 
simulations  which  rebuild  the  geometry  of  the  mechanical 
compression  by  solving  some  differential  equations  (e.g. 
compression  of  spherical  aggregates).  In  contrast,  the  second  class 
represents  a  purely  geometrical  approach  which  does  not  include 
any  physics  and  only  aims  to  realistically  rebuild  compressed 


microstructures.  Since  the  compression  of  microstructures  needs 
extensive  computational  time  for  physical  models  and  thus  limits 
the  volumes  or  resolution  that  can  be  handled,  we  focus  on  the 
second  class  of  compression  models  in  the  following.  Furthermore, 
within  this  class,  we  restrict  ourselves  to  uniaxial  compression.  In 
literature,  there  exist  some  models  for  uniaxial  compression  of 
fiber-based  materials,  cf.  [18-20  .  In  Ref.  [19],  a  compression  model 
is  proposed  which  mimics  compression  by  a  flat  stamp  only  for 
multi-layer  structures,  i.e.,  consisting  of  stacked  layers.  In  Refs. 
[18,20],  a  further  compression  model  is  presented  which  allows  to 
uniaxially  compress  a  given  3D  microstructure  of  some  fiber-based 
materials  by  a  flat  stamp  where  the  microstructure  is  represented 
by  a  3D  image.  So  far,  to  the  best  of  the  authors’  knowledge,  there 
exists  no  purely  geometrically  motivated  compression  model 
which  describes  uniaxial  compression  for  arbitrary  stamp  geome¬ 
tries.  However,  in  reality  the  GDL  is  compressed  by  a  flow-field 
stamp.  A  stamp  which  mimics  this  structure  is  shown  in  Fig.  1 
(right).  Note  that  there  are  dramatic  structural  differences  when 
compressing  the  GDL  either  with  a  flat  stamp  or  with  a  stamp 
possessing  a  flow-field  structure.  Thus,  there  is  a  need  for  more 
general  compression  models  in  order  to  analyze  the  influence  of 
compression  by  flow-field  stamps  on  transportation  properties  of 
the  GDL.  To  fill  this  gap,  we  present  a  model  which  is  able  to 
describe  all  kinds  of  uniaxial  compression,  in  particular  compres¬ 
sion  by  a  flow-field-like  stamp  as  occurring  for  GDL  in  operating 
PEMFC.  More  precisely,  given  the  input  of  a  3D  microstructure  of 
some  fiber-system,  the  model  uniaxially  compresses  the  system  of 
fibers  for  arbitrary  stamp  geometries  and  for  arbitrary  compression 
rates.  The  compression  model  is  constructed  by  means  of  a  two- 
stage  approach.  The  basic  idea  of  the  model  is  to  translate  the  fi¬ 
bers  in  the  direction  of  compression  according  to  a  vector  field  Vc 
which  depends  on  the  compression  rate  c  and  on  the  location  of  the 
fibers  within  the  material.  Subsequently,  an  avoidance  algorithm  is 
applied  to  the  translated  fiber  system  which  prevents  the  mutual 
penetration  of  fibers  and  preserves  their  lengths.  Thus,  the  pro¬ 
posed  compression  model  generalizes  the  model  described  in  Ref. 
[20]  in  two  ways:  1)  to  arbitrary  stamp  geometries,  and  2)  it  keeps 
non-overlapping  fibers  free  of  overlaps.  Moreover,  an  estimation 
technique  for  the  optimal  vector  field  in  order  to  fit  it  to  experi¬ 
mental  data  is  introduced  which  is  based  on  simulated  annealing. 

To  validate  the  model,  it  is  applied  to  the  3D  microstructure  of 
non-woven  GDL  in  PEMFC  gained  by  synchrotron  tomography  for 
four  different  compression  rates  ce  {0, 1 0, 20, 30}  (in  %).  Thereby,  we 
distinguish  between  two  kinds  of  compression:  1 )  using  a  flat  stamp, 
and  2)  a  stamp  possessing  a  flow-field  structure.  Note  that  when  the 
compression  model  is  applied  to  systems  of  non-overlapping  fibers, 
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Fig.  2.  Scheme  of  the  principal  experimental  set  up  at  the  tomographic  instrument. 

single  fibers  need  to  be  identified  in  order  to  separate  them  from  each 
other.  Such  a  single-fiber  extraction  can  be  done  efficiently  by  the 
algorithm  introduced  in  Ref.  [21  ].  The  compression  model  is  vali¬ 
dated  by  comparing  structural  characteristics  computed  for  exper¬ 
imentally  and  virtually  compressed  microstructures.  For  both  types 
of  stamps,  a  good  agreement  between  experimental  and  virtual 
compression  is  found. 

In  a  second  step,  we  apply  the  compression  model  to  virtual 
microstructures  of  uncompressed  non-woven  GDL,  which  are 
generated  according  to  a  stochastic  3D  simulation  model  intro¬ 
duced  in  Ref.  [22].  The  combination  of  the  stochastic  microstructure 
model  with  the  compression  model  allows  to  realistically  simulate 
compressed  GDL  on  an  arbitrary  large  scale  for  arbitrary 
compression  rates.  Since  the  access  to  tomography  images  is 
limited  or  often,  only  images  of  uncompressed  media  are  available, 
the  combination  of  both  models  enables  us  to  efficiently  generate  a 
large  number  of  virtual  (but  realistic)  microstructure  scenarios  for 
compressed  GDL.  When  employed  in  fuel  cell  simulation,  this 
model  provides  great  flexibility:  the  microstructure  of  GDL  can  be 
easily  adapted  to  different  degrees  of  compression  in  order  to  study 
the  influence  of  clamping  pressure  on  fuel  cell  operation.  Further¬ 
more,  it  allows  for  studying  different  flow-field  geometries  while 
the  volume  of  the  generated  GDL-structure  can  be  conveniently 
adjusted  to  meet  the  desired  requirements.  By  systematically 
varying  the  parameters  of  the  compression  model  (describing  the 


3D  microstructure  of  compressed  GDL),  new  virtual  GDL  can  be 
generated  and  in  combination  with  computational  transport  sim¬ 
ulations,  the  microstructure-functionality  relationship  of  the  GDL 
in  operating  PEMFC  can  be  investigated.  In  this  way,  the  systematic 
design  of  morphologies  with  improved  functionality  based  on 
computer  experiments  can  be  achieved. 

The  paper  is  organized  as  follows.  Section  2  briefly  describes  the 
material  investigated  in  this  paper  and  the  set-up  for  the  genera¬ 
tion  of  3D  image  data.  In  Section  3,  the  compression  model  is 
introduced  and  the  estimation  technique  for  its  parameter  is  pre¬ 
sented.  Subsequently,  in  Section  4,  the  compression  model  is 
applied  to  synchrotron  data  of  non-woven  GDL  for  different 
compression  rates,  where  we  distinguish  between  two  kinds  of 
compression,  i.e.,  using  a  flat  stamp  and  a  stamp  possessing  a  flow- 
field  structure.  Section  5  combines  a  stochastic  simulation  model 
for  uncompressed  non-woven  GDL  with  the  compression  model 
introduced  in  the  present  paper.  Finally,  Section  6  summarizes  the 
results  and  provides  a  short  outlook  regarding  possible  future 
research. 

2.  Imaging  and  material 

A  prominent  example  of  a  fiber-based  material  is  the  GDL  in 
PEMFC  where  the  consideration  of  compressed  GDL  is  of  great 
importance  since  the  GDL  is  compressed  in  the  operating  state  of 
PEMFC.  Therefore,  3D  image  data  of  uncompressed  and  compressed 
non-woven  GDL  are  needed  23],  where  we  distinguish  between 
compression  by  a  flat  and  by  a  flow-field-like  stamp,  see  Figs.  3  and 
4.  Thus,  in  this  section,  we  briefly  discuss  the  imaging  technique 
used  in  Ref.  [23]  to  record  compressed  microstructures  and  we  give 
some  details  about  the  material  considered.  Finally,  we  shortly 
discuss  an  algorithm  which  automatically  extracts  a  system  of 
single  fibers  from  tomographic  3D  images.  This  algorithm  allows 
the  direct  application  of  the  compression  model,  introduced  in 
Section  3,  to  experimental  image  data  of  fiber-based  materials. 

2 A.  Imaging  technique 

The  imaging  experiments  were  performed  at  the  Helmholtz- 
Zentrum  Berlin  using  the  synchrotron  tomography  station  of  the 
BAMline  which  is  located  at  the  synchrotron  source  BESSY  (Berlin / 
Germany)  [24].  The  high  brilliance  of  the  synchrotron  X-ray  beam 
allows  producing  tomographic  data  of  superior  quality  in  terms  of 


Fig.  3.  3D  image  data  gained  by  synchrotron  tomography  of  non-woven  GDL  compressed  by  a  flat  stamp  with  compression  rates  c  =  0%  (top  left),  c  =  10%  (top  right),  c  =  20% 
(bottom  left)  and  c  =  30%  (bottom  right). 
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Fig.  4.  3D  image  data  gained  by  synchrotron  tomography  of  non-woven  GDL  compressed  by  a  flow-field  stamp  with  compression  rates  c  =  0%  (top  left),  c  =  10%  (top  right),  c  =  20% 
(bottom  left)  and  c  =  30%  (bottom  right).  The  red  (yellow)  marked  fibers  are  located  underneath  the  channel  (the  rips  of  the  flow-field  stamp).  (For  interpretation  of  the  references 
to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


resolution  and  contrast  [23,25,26].  This  way,  smallest  details  of  the 
three-dimensional  GDL  structure  such  as  shape  and  course  of  in¬ 
dividual  carbon  fibers  can  be  accurately  rendered.  The  principal  set 
up  of  the  tomographic  measurement  is  displayed  in  Fig.  2. 

A  W-Si  double  multilayer  monochromator  with  an  energy 
resolution  of  A  E/E  =  10  2  was  used  to  obtain  a  monochromatic  X- 
ray  beam.  The  beam  energy  was  adjusted  to  15  keV  in  order  to 
achieve  optimal  contrast  for  fibers.  Images  were  captured  with  a 
4008  x  2672  pixel2  camera  set-up  (PCO  camera  with  a  Gadox 
scintillator  screen)  rendering  a  field  of  view  of  3.3  x  2.2  mm2  with  a 
pixel  size  of  0.876  pm  and  a  respective  physical  spatial  resolution 
of  about  2  pm  [27].  Circular  GDL  samples  with  diameters  of  3  mm 
were  adjusted  in  the  sample  holder  and  mounted  on  a  translation/ 
rotation  unit.  For  each  tomogram,  samples  were  rotated  in  equi¬ 
distant  steps  covering  an  angular  range  of  180°.  A  radiographic  set 
of  1500  projections  and  500  flatfields  was  taken  and,  subsequently, 
reconstructed  to  a  3D  volume.  The  exposure  time  is  2.5  s  plus  1.7  s 
read-out  for  a  single  radiograph  adding  up  to  an  acquisition  time  of 
140  min  for  the  complete  tomogram.  The  measurements  were 
performed  using  a  dedicated  compression  device  as  sample  holder 
which  provides  well-defined  compression  conditions  and  meets 
the  instrumental  requirements,  e.g.  in  terms  of  sufficient  beam 
transmission.  The  compression  device  is  able  to  compress  GDL  by 
different  stamp  types  (see  Figs.  3  and  4),  in  particular,  for  a  flat 
stamp  and  a  stamp  consisting  of  a  flow-field  structure,  where  the 
width  of  the  flow-field  canal  is  790  pm. 

2.2.  Description  of  material 

The  application  of  the  compression  model,  introduced  in  the 
present  paper,  is  based  on  tomographic  data  of  the  GDL  type  FI2315, 
a  non-woven,  carbon  fiber-based  material  produced  by  the  com¬ 
pany  Freudenberg  FFCCT.  This  material  does  not  contain  binder,  wet 
proofing  agents  nor  a  micro  porous  layer  (MPL).  FI2315  is  used  as 
standard  material  for  many  fuel  cell  applications.  Note,  however, 
that  it  is  usually  equipped  with  wet  proofing  agent  and/or  MPL  in 
order  to  adapt  its  physical  properties  to  the  requirements  of  the 
respective  application. 

2.3.  Fiber  extraction  algorithm 

The  experimental  3D  image  data  gained  by  synchrotron  to¬ 
mography  displays  the  (uncompressed  and  compressed)  non- 


woven  GDL  by  a  matrix  of  voxels.  In  the  following,  we  extract  the 
systems  of  fibers  from  the  tomographic  3D  image  data  of  uncom¬ 
pressed  GDL.  Thereby,  the  complexity  of  data  is  reduced  since  it  is 
represented  by  thousands  of  fibers  (approximated  by  polygonal 
tracks  which  are  given  in  Euclidean  coordinates)  instead  of  millions 
of  voxels.  Furthermore,  the  input  data  of  the  compression  model 
introduced  in  this  paper  are  uncompressed  fiber-systems  and  thus 
it  can  be  directly  applied  to  the  extracted  fibers  of  uncompressed 
non-woven  GDL. 

We  briefly  describe  an  algorithm  to  automatically  detect  single 
fibers  from  3D  tomographic  data  of  fiber-based  materials  which  has 
been  proposed  in  Ref.  [21].  The  algorithm  yields  the  extracted 
center  lines  of  individual  fibers  which  are  approximated  by 
polygonal  tracks.  Since  in  experimental  data  there  occur  irregu¬ 
larities  like  noise  or  binarization  artefacts  it  is  only  possible  to 
extract  parts  of  the  center  lines.  Then,  a  stochastic  optimization 
algorithm  connects  these  fragments  of  the  center  lines  to  each 
other  such  that  the  complete  fibers  are  correctly  reconstructed.  For 
further  information  see  Ref.  [21].  Exemplarily,  we  have  applied  the 
algorithm  to  experimental  image  data  of  uncompressed  non- 
woven  GDL.  The  experimental  data  and  its  extracted  fiber  system 
are  in  an  excellent  optical  accordance,  see  Fig.  5. 

3.  Compression  model  for  fiber-based  materials 

In  this  section,  we  introduce  a  parametrized  model  that  de¬ 
scribes  the  3D  microstructure  of  compressed  materials.  Given  the 
input  of  an  uncompressed  3D  microstructure  of  a  two-phase  ma¬ 
terial  consisting  of  an  uncompressed  solid  phase  and  a  pore  phase, 
the  model  compresses  the  solid  phase  in  a  uniaxial  direction  for 
arbitrary  compression  rates.  When  speaking  of  a  3D  microstructure, 
we  refer  to  a  segmented  tomographic  3D  image  or  synthetic  3D 
images  drawn  e.g.  from  stochastic  simulation  models.  To  make  the 
model  more  accessible,  we  constrict  it  and  its  range  of  application 
to  materials  consisting  of  non-overlapping  systems  of  fibers.  But 
note  that  by  slight  modifications,  the  model  can  be  extended  to 
other  types  of  materials,  including  composite  materials.  The 
compression  model  introduced  in  the  following  should  possess 
several  properties.  In  particular,  the  model  should 

1)  describe  different  uniaxial  compression  processes  on  the 

microstructure  of  fiber-based  materials  realistically,  where  it 

should  describe  the  compression  of  fiber-based  materials 
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Fig.  5.  3D  image  of  uncompressed  GDL  (left)  and  its  extracted  system  of  single  fibers  (right). 


accomplished  by  different  stamp  geometries,  e.g.  by  a  flat  stamp 
or  by  a  stamp  with  a  flow-field  structure,  see  also  Fig.  1. 

2)  preserve  the  non-overlapping  of  fibers. 

3)  describe  compressed  microstructures  of  fiber-based  materials 
with  continuous  compression  rates  ce  [0,cmax),  where 

Cmax  <  1 

is  the  largest  compression  rate  such  that  it  is  possible  to  pre¬ 
serve  a  non-overlapping  fiber-system  after  applying  the 
compression  model. 

4)  be  a  volume  preserving  operation,  i.e.,  the  total  volume  of  fibers 
should  remain  constant  under  compression. 

In  the  following,  we  first  introduce  the  compression  model  for 
fiber-based  materials.  Then,  a  fitting  procedure  of  its  parameters  is 
discussed  when  the  uncompressed  and  several  compressed  mi¬ 
crostructures  of  a  fiber-based  material  are  given. 


3.1.  Model  description 

We  now  introduce  a  compression  model  which  describes  uni- 
axially  compressed  microstructures  of  fiber-based  materials. 
Throughout  this  paper,  we  assume  that  microstructures  are  com¬ 
pressed  in  a  uniaxial  direction,  say,  the  z-direction.  Consequently, 
we  assume  that  the  fibers  are  only  shifted  in  the  direction  of 
compression,  but  the  magnitude  of  the  shift  may  vary  in  space.  To 
motivate  the  assumption,  consider  Fig.  6  displaying  a  cross-section 
of  a  fiber  system  (GDL)  in  an  uncompressed  state  (left)  and  a 
compressed  state  (right).  It  can  clearly  be  seen  that  fibers  are 
mainly  shifted  in  z-direction.  Note  that  due  to  binarization  artefacts 
and  irregulations  of  the  imaging  of  the  material,  it  is  implausible  to 
obtain  a  perfect  translation  in  z-direction  even  if  100%  of  the  fibers 
were  shifted  in  z-direction.  On  the  other  hand,  it  is  intuitively  clear 
that  fibers  are  translated  to  a  large  part  only  in  z-direction  when 
considering  the  case  of  uniaxial  compression.  For  these  reasons,  we 
neglect  possible  translations  of  fibers  in  other  directions. 

The  compression  model  is  constructed  by  means  of  a  two-stage 
approach.  First,  the  fibers  are  translated  in  the  direction  of 


compression  according  to  a  vector  field  Vc  depending  on  the 
compression  rate  ce  [0,  cmax).  Note  that  each  fiber  is  represented  by 
a  polygonal  track.  To  the  starting  and  endpoints  of  each  segment  of 
the  polygonal  tracks,  the  vector  field  assigns  a  transition  vector, 
which  depends  on  the  location  of  the  point  within  the  material.  The 
translation  of  the  fibers  according  to  the  vector  field  Vc  may  yield 
overlaps  in  the  fiber  system  and  increase  or  decrease  the  length  of 
individual  fibers.  Therefore,  in  a  second  step,  we  apply  an  iterative 
avoidance  algorithm  [28]  to  the  translated  fiber  system  which 
transforms  the  translated  fiber  system  to  a  non-overlapping  fiber 
system,  where  the  length  and  curvature  of  the  fibers  are  preserved. 

3.1.1.  Translation  algorithm 

Let  S  =  Uj“1Fj©B(0,  r)c^  x  U  x  [0,zmax]  be  a  system  of  fibers 
(representing  a  fiber-based  material)  and  let  FJ,  i  =  1,  2,  ...  be 
polygonal  tracks  approximating  single  fibers,  i.e.,  FJ  is  a  system  of 
line  segments  defined  by 

rii 


where 

fiy  =  (P?),P?>«3X2 

is  the  jth  line  segment  of  the  ith  fiber  consisting  of  starting  and 
endpoints  p^,pj!W  with  Euclidean  length  \Hy\  =  r/4.  More¬ 
over,  zmax  >  0,  r  defines  the  constant  radius  of  the  fibers,  B(0,r)  the 
(three-dimensional)  sphere  with  origin  as  center  point  and  radius  r, 
and  ©  denotes  the  Minkowski  addition.  Thus,  we  identify  the 
microstructure  of  a  fiber-based  material  with  the  union  set  of  all 
fibers,  where  each  fiber  is  represented  by  a  system  of  connected 
line  segments.  Given  this  representation,  we  can  now  describe  the 
compression  of  fibers  in  more  detail.  Let  Vc :  U  x  U  x  [0,  zmax]  ->  U 
be  a  continuous  vector  field  with  the  constraints  that  Vc(x,y,0)  =  0 
for  all  (xj)tgR2,  and  Vc(x,y,z)  <  z  for  all  (x,y)Te  U2  and  0  <  z  < 


Fig.  6.  Cross-section  view  on  non-woven  GDL  in  uncompressed  state  (left)  and  in  compressed  state  with  c  =  20%. 
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Zmax*  In  the  first  stage,  the  system  of  fibers  E  is  translated  according 
to  the  vector  field  Vc  which  yields  the  translated  system  of  fibers, 
denoted  by  S',  where 

S'  =  U  F^©E(0,r),  with  F[  =  U  and 

Thus,  the  translation  of  the  fiber  system  is  described  by  trans¬ 
lating  the  starting  and  endpoints  of  each  line  segment  of  the  fibers. 
In  a  second  step,  an  avoidance-algorithm  is  applied  to  the  trans¬ 
lated  system  of  fibers  S',  in  order  to  de-overlap  the  system  of  fibers. 
The  resulting  microstructure  is  denoted  by  S".  Note  that  the 
compressed  microstructure  S"(c)  as  well  as  the  translated 
fiber  system  S'(c)  are  located  in  the  set  Wc,  where 
Wc  =  {(x,y,z)TeR3  :  (x,y)TeR2,ze[0,zmax  -  yc(x,y,zmax)]}.  In 
the  following,  we  shortly  discuss  the  iterative  avoidance  algorithm 
introduced  in  Ref.  [28]. 

3.1.2.  Iterative  avoidance  algorithm 

To  obtain  a  system  of  compressed  and  non-overlapping  fibers, 
the  translated  fiber  system  S'  described  in  Section  3.1.1  is  further 
transformed.  Recall  that  r  denotes  the  thickness  of  the  fibers.  The 
principle  idea  is  to  first  represent  each  fiber  by  a  chain  of  (over¬ 
lapping)  spheres  where  the  spheres  have  the  same  radius  r  equal  to 
the  fiber  thickness  and  their  midpoints  are  located  equidistantly 
(with  distances  r/4  from  each  other)  on  the  line  segments  of  the 
fibers,  see  also  Fig.  7  (left  and  center).  Using  this  representation,  the 
basic  idea  is  to  translate  the  midpoints  of  the  spheres  such  that  no 
spheres  from  different  fibers  overlap,  see  Fig.  7  (right). 

In  more  detail,  in  every  iteration  step  i,  each  midpoint  M^e  U3 
of  those  spheres  is  translated  by  a  vector  t  =  t\  +  t2 ,  i.e., 

a4'+1)  =  M®  +  t  =  M®  +h  +  t2. 

The  translation  function  t\  separates  overlapping  spheres  which 
belong  to  different  chains  of  spheres  by  moving  them  away  from 
each  other,  and  the  second  translation  t2  preserves  the  distance  and 
the  angle  between  neighboring  spheres  belonging  to  the  same 
chain  of  spheres  (and  thus  the  length  and  curvature  of  the  fibers). 
For  the  application  of  the  compression  model  within  the  present 
paper,  the  function  t\  is  weighted  stronger  than  t2  as  proposed  in 
Ref.  [28],  i.e.,  to  keep  fibers  non-overlapping  is  more  important 
than  preserving  their  curvatures  and  lengths.  Note  that  if  consid¬ 
ering  compression  of  fiber-based  materials  where  the  elasticity  of 
fibers  is  strongly  restricted,  i.e.,  large  bending  of  fibers  is  not  real¬ 
istic  without  cracking,  the  weighting  of  functions  t\  and  t2  can  be 
adopted  accordingly.  In  this  case,  the  fibers  would  not  be  allowed  to 
exceed  a  certain  bending.  If  after  an  iteration  step  i,  a  midpoint 
=  (mfy,  mfy,  m£z)T  is  not  located  in  Wc,  it  is  set  to  the  point 
given  by  the  closest  distance  with  respect  to  the  z-coordinate,  i.e.,  if 
M®£WC  we  set  M®  =  argmin(mm  /)eWJz  -  m®z|.  For  further 


information  regarding  the  iterative  avoidance  algorithm,  the  reader 
is  referred  to  [28].  The  final  system  of  non-overlapping  fibers  is 
denoted  by  E"(c). 

3.2.  Fitting  of  vector  field  Vcfor  fixed  compression  rate  c 

So  far,  the  compression  model  has  been  introduced  which  is 
uniquely  determined  by  the  parameters  of  the  continuous  vector 
field  Vc.  It  remains,  however,  unclear  which  choice  of  Vc  is  most 
suitable  to  describe  compression  of  fiber-based  materials  in  a 
realistic  way.  Therefore,  in  this  section,  we  introduce  a  fitting  pro¬ 
cedure  which  estimates  the  continuous  vector  field  Vc  of  the 
compression  model  described  in  Section  3.1  for  any  given  value  of 
ce  [0,  cmax).  For  this  procedure,  3D  image  data  of  one  and  the  same 
fiber  system  in  both,  an  uncompressed  and  a  compressed  state,  is 
required.  The  estimation  technique  is  based  on  the  fact  that  the 
total  volume  of  the  fiber  system  must  remain  constant  under 
compression,  since  fibers  cannot  overlap.  More  precisely,  the  idea  is 
to  choose  a  vector  field  such  that  the  translated  fiber  system  E" 
matches  the  compressed  fiber  system  observed  in  experimental 
data  as  good  as  possible.  Since  we  cannot  expect  a  perfect  match  of 
these  fiber  systems  due  to  e.g.  binarization  artefacts,  this  procedure 
optimizes  the  match  of  local  distributions  of  fibers  in  the  virtually 
compressed  and  the  experimentally  compressed  microstructures. 
In  particular,  we  subdivide  the  uncompressed  and  compressed 
microstructures  in  small  (pairwise  disjoint)  cuboids  and  minimize 
the  absolute  values  of  the  differences  of  the  total  fiber  volumes  in 
corresponding  cuboids  for  uncompressed  and  compressed  data.  To 
solve  this  optimization  problem  we  use  simulated  annealing,  a 
well-known  stochastic  optimization  algorithm. 

3.2.1.  Optimization  problem 

Let  us  assume  that  we  consider  the  uncompressed  microstruc¬ 
ture  in  the  3D  volume  W0  SR3  and  the  compressed  microstructure 
in  W\  [R3.  Furthermore,  we  assume  that  the  dimension  in  x-  and 
y-direction  of  W0  and  W\  are  the  same  but  by  the  uniaxial 
compression  the  dimensions  in  z-direction  are  allowed  to  vary,  see 
Fig.  8a.  In  a  first  step,  the  observation  windows  Wo  and  W\  are 
subdivided  in  the  same  number  of  small  disjoint  cuboids,  where  the 
side  length  of  the  cuboids  are  equal  in  x-  and  y-direction  for  Wo  and 
W\  but  differ  in  z-direction,  see  Fig.  8b.  Note  that  the  decomposi¬ 
tion  of  Wo  and  W\  is  much  finer  than  displayed  in  Fig.  8b. 

Given  the  decomposition  of  Wo  and  W\,  we  now  consider  a 
single  pair  of  stacks  of  cuboids,  as  displayed  in  Fig.  8c  (left).  We 
denote  both  stacks  by  S'  =  (Zg,...,z„)  and  S  =  (zo,  ...»  zn)  with 
corresponding  total  volumes  of  fibers  within  the  cuboids  given  by 
(V[, V')  and  (Vi(z0,  ...,  z„),  ...,  Vn{z0,  z„)),  respectively.  Since 
we  assume  that  fibers  are  only  translated  in  z-direction,  the  total 
fiber  volume  of  both  stacks  must  be  the  same.  Thus,  the  fibers  have 
to  be  compressed  in  such  a  way  that  the  fiber  volume  of  each 
cuboid  in  the  uncompressed  microstructure  is  the  same  as  the  fiber 
volume  of  the  corresponding  cuboid  in  the  compressed 


Fig.  7.  Schematic  2D  visualization  of  a  fiber  system  (left)  being  represented  by  overlapping  spheres  (center)  and  subsequently  translated  to  a  system  without  overlaps  between 
different  fibers  (right). 
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a)  Cross-sectional  view  on  fiber  system 


b)  Decomposition  of  observation  windows  W0  and  Wx  in  cuboids 


uncompressed 


compressed 


C)  For  each  stack  independently:  Optimize  decomposition  of  d)  Result: 

cuboids  such  that  volumes  of  cuboids  coincide 


Resulting  vector  field: 
length  of  vectors  are 
given  by  z*j-z'j 

Vector  field  is  further 
processed  to  continous 
field  Vc 


Fig.  8.  Fitting  the  vector  field  of  the  compression  model  for  fixed  c. 


microstructure.  In  other  words,  the  ‘optimal’  decomposition  S* 
should  fulfill  the  equality 

S*  =  arg  min  V \Vt{Zo,  -  V{\2.  (1) 

0  =  z0<...<z„  J 

Thus,  we  have  formulated  an  optimization  problem  and  explain 
later  on  (in  Section  3.2.2),  how  to  solve  it.  Note  that,  since  we  as¬ 
sume  that  fibers  are  only  translated  in  z-direction,  it  is  sufficient  to 
optimize  each  vertical  stack  of  cuboids  separately,  see  Fig.  8c.  For 
the  moment,  let  us  assume  that  we  have  already  obtained  the 
optimal  decomposition  S*  =  (zg,  ...5z*).  Then  this  gives  us  detailed 
information  about  the  local  compression:  Each  starting  and  end¬ 
points  of  the  fiber  segments  (forming  polygonal  tracks)  of  the  un¬ 
compressed  fiber  system  that  is  located  at  point  ( Xj,yk,z ?)  should 
be  translated  by  the  vector  given  by  (0, 0,z-  -  z- ) T  according  to  the 
optimal  decomposition  5*.  Thus,  the  result  of  the  optimal  decom¬ 
position  S*  is  a  set  of  translation  vectors  given  on  discrete  points  in 
the  observation  window  Wo,  see  Fig.  8d. 

To  compute  the  ’optimal’  vector  field  Vc,  we  fit  a  continuous  and 
parametric  function  to  these  translation  vectors  by  non-linear 
regression  [29].  In  particular,  we  first  consider  a  suitable  class 
of  parametrized  continuous  vector  fields  {Vc(x,y,z)  : 
Vc((x,y,  0)  =  0,  and  Vc(x,y,z)  <z  for  all  (x,y)Te[R2}  such  that  no 
transitions  beyond  the  observation  window  are  made.  Then,  the 
parameters  of  the  vector  field  are  estimated  by  non-linear  regres¬ 
sion  methods.  For  this,  we  use  the  surface  fitting  tool  implemented 
in  Matlab  [30].  To  check  the  quality  of  how  the  estimated  contin¬ 
uous  vector  field  fits  the  translation  vectors  given  on  discrete  points 
in  space,  we  consider  the  coefficient  of  determination  R2e  [0, 1).  It  is 
a  measure  for  goodness-of-fit  in  terms  of  non-linear  regression.  The 
closer  R2  is  to  1  the  better  the  regression  model  describes  the  given 
data.  Note  that  R2  >  0.9  represents  a  very  good  fit  [29],  where  a  very 
good  data  fit  means  that  the  local  distribution  of  fibers  in  the 
virtually  compressed  structure  coincides  with  the  distribution 
found  in  the  image  data  of  the  experimentally  compressed 
structure. 

3.2.2.  Simulated  annealing 

In  order  to  solve  the  optimization  problem  described  in  (Section 
1 ),  we  use  simulated  annealing,  a  well-know  stochastic  optimization 


algorithm  which  numerically  solves  complex  global  optimization 
problems  by  means  of  Markov-Chain-Monte-Carlo  simulation  [31]. 
Roughly  speaking,  the  basic  idea  of  the  simulated  annealing  algo¬ 
rithm  is  to  start  with  some  initial  configuration,  say  (z^0) ,  ...,z^) 
and  then  randomly  change  one  component  of  the  vector.  The 
change  is  accepted  if  the  cost  function  (which  measures  the  devi¬ 
ation  of  the  fiber  volumes  within  the  uncompressed  and  com¬ 
pressed  cuboids  and  which  we  aim  to  minimize)  gets  smaller. 
Otherwise,  the  change  is  only  accepted  with  a  certain  probability 
and  the  next  iteration  step  is  started.  The  interesting  aspect  of 
simulated  annealing  is  that  it  allows  transitions  which  increase  the 
cost  function.  This  makes  an  ‘escape’  from  local  minima  possible. 
The  probability  for  accepting  transitions  that  increase  the  cost 
function  decreases  with  an  increasing  number  of  steps  and  is 
governed  by  a  so-called  cooling  schedule. 

We  now  explain  the  simulated  annealing  algorithm  in  more 
detail.  Let  the  side  length  in  z-direction  of  the  cuboids  of  Wo  be 
given  by  Zq  =  0  <  z\  <  ...  <z'n  =  zmax  and  the  corresponding  side 
length  of  the  fixed  decomposition  of  W\  be  given  by 
Zo  =  0  <  Z\<  ...  <  zn.  Furthermore,  we  denote  by  MSE(zo, zn)  the 
mean  squared  deviation  of  the  volumes  computed  for  corre¬ 
sponding  cuboids  given  the  decomposition  of  Wo  by  (zg,  . . . ,  zg)  and 
Wi  by  (z0,  ...,z„),  i.e., 


MSE(z0,...,z„)  =  ^|Vi(z0,...,z„)-  V!\2. 

1=1 

The  simulated  annealing  algorithm  has  three  parameters 
T  >  0,  Mg  N  and  g<=  (0,1).  Note  that  the  tripel  (T,n,g)  is  called  cooling 
schedule,  see  Ref.  [31  ].  The  algorithm  works  as  follows,  where  we 
putz<,0)  =  z0,...,z*,0)  =  Zn- 


(a)  Set  m  =  1  and  repeat  the  steps  (b)  to  (d)  until  m  =  M  then  go 
to  (e). 


(b)  Choose  randomly  one  z\m  ^e{z^m 


>-i) 


m 


=  zjm  ^  for  all  je{l,  ...,n}/{z}. 
(c)  Modify  z|m_1^  by  adding 


}  and  put 


random 


X~LT(— 0.01,0.01), 

Jjn- 1)  ,  v  ^  Jjn- 1) 


i.e., 


zm  =  z(m- 1)  +Xm 


number 
If 


z;-  ])  +X  >  or  zr-'>  +X  <  zr~l),  zf "  =  zr-‘>  and 

repeat  step  (c). 
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(d)  If  MSE(z^m\...,z^m))  <  MSE(z£"“1),...,z<1m“1))  accept  the 
modification,  otherwise  accept  only  with  a  certain  proba¬ 
bility  which  is  given  by 

exp(-  (MSE(z<,m),..„z<m))  -MSE(z<,m-1),...,zf-1)))/T). 

Set  m  =  m  +  1  and  continue  with  (b). 

(e)  Set  T  =  g  T.  If 


|MSE(z^...,zr)-MSE(zr,...,zr)| 

_ ! _ X-i _ L _ X-i _ IA _ <  o  0001 

max{  MSE  (z<,0) , . . . ,  z<0) ) ,  MSE  (z<,M) , . . . ,  z(nM) )  } 

the  algorithm  terminates,  otherwise  set  z ^  =  z^ , . . . ,  zj,0)  =  z^M) 
and  go  to  (a). 

As  initial  value  for  T  we  choose  T  =  1  10~6,  which  is  a  value 
where  enough  changes  are  accepted.  The  number  M  of  iterations 
per  step  is  chosen  equal  to  M  =  10,000,  and  the  factor  g  is  put  equal 
to  g  =  0.98  as  recommended  in  Ref.  [31  ]. 

3.3.  Validation  of  fitting  procedure 

In  this  section,  we  analyze  the  goodness  of  the  fitting  procedure 
for  the  continuous  vector  field  Vc  for  fixed  ce  [0,  cmax)  introduced  in 
Section  3.2.  Therefore,  we  generate  uncompressed  synthetic  3D  fiber 
systems,  using  the  stochastic  model  presented  in  Ref.  [22].  Then,  for 
a  given  vector  field  Vc ,  we  virtually  compress  the  synthetic  fiber 
system  using  the  compression  model  introduced  in  Section  3.1  and 
subsequently  estimate  the  vector  field  Vc  from  these  uncompressed 


virtual  microstructures,  as  described  in  Section  3.2.  Then,  we 
compare  the  given  vector  field  Vc  to  its  estimated  counterpart  Vc. 

The  parameters  of  the  stochastic  fiber  model  introduced  in  Ref. 
[22]  are  chosen  such  that  the  virtually  generated  fiber  systems 
resemble  the  fiber  systems  of  experimental  (uncompressed)  non- 
woven  GDL,  see  Section  2.  Note  that  the  uncompressed  (virtual) 
microstructures  are  simulated  on  the  observation  window 
W0  =  [0,800]  x  [0,800]  x  [0,200]  (i.e.,  dimensions  of  generated  GDL 
structures),  see  also  Fig.  9  (left). 

The  compressed  microstructures  are  obtained  by  applying  the 
compression  model  presented  in  Section  3.1,  for  a  given  vector  field 
Vc  to  the  simulated  (uncompressed)  fiber  system.  In  particular,  we 
consider  two  different  vector  fields  and  vj?\  where  the  first 
one  is  chosen  as  a  linear  vector  field  V^\x,y,z)  =  0.5  z  which  only 
depends  on  the  z-component,  whereas  the  second  one  is  a  vector 
field  V^2\x1y,z)  =  z- (0.3  -  |x  -  400|  0.0015)  depending  on  the  x- 
and  z-components.  Subsequently,  we  apply  the  compression  model 
for  these  two  vector  fields  to  the  synthetic  uncompressed  data  and 
compare  the  fitted  vector  fields  to  the  actual  ones.  Note  that  both 

vector  fields  and  fulfill  the  properties  of  a  vector  field  for 
the  compression  model,  i.e.,  both  are  continuous,  and  it  holds  that 
V®(x,y,  0)  =  0  and  V^(x,y,z)  <  z  for  ie{l,2}.  The  uncompressed 
microstructure  generated  by  the  stochastic  3D  simulation  model  as 
well  as  the  compressed  microstructures  gained  by  applying  the 

compression  model  for  the  vector  fields  and  V ^  to  the  un¬ 
compressed  microstructures  are  visualized  in  Fig.  9. 

Recall  that  by  means  of  the  fitting  technique  proposed  in  Section 
3.2,  estimates  of  the  vector  field  V$p  at  discrete  points  in  the 
observation  window  W0  are  computed.  Using  the  surface  fitting 
tool  implemented  in  Matlab  [30],  a  given  non-linear  surface  can  be 
fitted  to  these  estimates. 


Fig.  9.  Left  column:  realizations  of  two  uncompressed  microstructures  drawn  from  the  stochastic  3D  simulation  model;  right  column:  compression  model  applied  to  the  un¬ 
compressed  microstructures  (left)  for  the  vector  fields  (top)  and  (bottom). 
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Fig.  10.  The  fitted  vector  fields  (left)  and  V (right)  displayed  by  the  surface,  where  the  blue  points  are  estimates  of  the  vector  field  V®  at  discrete  points  in  the  observation 
window  W0  obtained  by  the  simulated  annealing  algorithm.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this 
article.) 


In  Fig.  10,  the  result  of  the  fitting  procedure  can  be  seen.  For  the 
first  vector  field  V^\  we  get  an  estimated  vector  field  of 

V^\x:y:z)  =  0.505  z  and  for  the  second  vector  field  an  esti- 

—  (2) 

mationofVc  (x,y,z)  =  z- (0.321  -  \x  -  400|  0.00169)  is  obtained. 

i  ~(i) 

Thereby,  the  coefficients  of  determination  Rz  =  0.977  for  Vc  and 

?  -(2) 

R  =  0.984  for  Vc  are  achieved.  Thus,  by  the  fitting  procedure 
introduced  in  Section  3.2,  it  is  possible  to  estimate  the  vector  field 
of  a  uniaxial  compression  almost  perfectly. 

3.4.  Arbitrary  compression  rates  c 

In  Sections  3.1  and  3.2,  a  toolbox  has  been  introduced  which 
consists  of  a  model  to  describe  microstructures  of  fiber-based 
materials  under  compression  and  a  procedure  to  estimate  the 
vector  field  Vc  of  the  compression  model  for  a  specified  compres¬ 
sion  rate  c.  Therefore,  the  uncompressed  and  one  single  com¬ 
pressed  microstructure  of  a  fiber-based  material  have  to  be  given. 
Thus,  it  is  possible  to  describe  the  microstructural  changes  from  the 
uncompressed  to  the  compressed  microstructure  for  one  single  rate 
of  compression  c.  Note  however  that  at  the  present  stage  it  is  not 
yet  possible  to  describe  the  compression  process  for  arbitrary  rates 
of  compression  ce  [0,  cmax). 

The  basic  idea  to  obtain  a  compression  model  for  arbitrary  rates 
of  compression  is  as  follows:  First,  let  us  suppose  that  an  uncom¬ 
pressed  and  m  compressed  microstructures  with  compression  rates 
ce  {c\ , . ..,  cm}  are  available.  For  each  ce  {c^ , ...,  cm},  we  estimate  the 
corresponding  vector  field  Vc  by  means  of  the  fitting  technique 
which  was  introduced  in  Section  3.2.  Assuming  that 
{’ Vc ,  ce  [0,  cmax)}  is  a  function  of  c,  this  gives  us  m  points  of  support 
of  Vc  (i.e.  {VCl, ...,  VCm })  and  by  subsequently  applying  regression 
techniques  we  obtain  a  set  of  vector  fields  (Vc,  ce  [0,  cmax)}. 

More  precisely,  the  compression  of  the  GDL  is  managed  by  the 
family  of  continuous  vector  fields  {Vc,  ce  [0,  cmax)}  which  is 

described  by  some  parameters  (pi,  ...,pn),  i.e.,  Vc  =  v^Pl,'",Pn\  Since 
the  vector  fields  depend  on  the  compression  rate,  we  interpret  the 


Table  1 

Estimated  coefficients  of  p,(c)  =  a,-(exp(bjc)  -  exp(d,c))  for  ie{2, 6}. 


i  =  2 

i  =  3 

i  =  4 

i  =  5 

i  =  6 

dj 

-0.235 

0.041 

0.784 

-1.282- 10-5 

-0.305 

bi 

-1.911 

0.347 

-7.583 

-1.150 

-1.666 

di 

-16.14 

0.316 

-7.627 

-20.15 

-1.666 

parameter  vector  (pi, . . .,  pn)  as  a  function  of  the  compression  rate  c, 
i.e. 

(Pl>-,Pn)T  :  [0,1]  —  Rn  (2) 

c«(pi(c),...,p„(c))T  1  J 

For  each  pair  of  uncompressed  and  compressed  microstructure, 
the  parameters  (pi,  ...,  pn)  of  the  same  parametric  vector  field 
v(Pi,...,pn)  are  fated  using  the  estimation  method  introduced  in 
Section  3.2.  This  yields  m  +  1  points  of  support  for  the  function  (pi, 
...,pn),  i.e.,  we  obtain  (Pi(c),  ...,pn(c))  force{0,Ci,  ...,cm}.  Based  on 
these  fitted  values  we  describe  the  dependency  of  the  compression 
rate  c  on  the  individual  parameters  p2  by  a  function  p2(c),  where  p2(c) 
is  computed  by  means  of  non-linear  regression.  The  vector  field  Vc 
which  uniquely  determines  the  compression  model  is  then  given 

by  v<Pi(c)’-«P»(c5). 

4.  Application  to  experimental  data 

To  begin  with,  we  apply  the  compression  model  introduced  in 
Section  3  to  the  experimental  image  data  of  non-woven  GDL  in 
PEMFC  which  have  been  described  in  Section  2.  In  particular,  as 
data  basis,  we  consider  experimental  3D  image  data  of  non-woven 
GDL  without  compression  (c  =  0)  and  for  three  different 
compression  rates  ce  {10, 20, 30}  (in  %).  The  goal  in  this  section  is  to 
use  the  model  introduced  in  Section  3  in  order  to  describe  the 
compression  of  fiber  systems  of  non-woven  GDL  in  PEMFC  for 
arbitrary  rates  of  compression  ce  [0,  cmax).  Finally,  we  validate  the 
compression  model  for  arbitrary  rates  of  compression  by 
comparing  structural  characteristics  computed  for  experimentally 
and  virtually  compressed  microstructures. 

4.1.  Estimation  of  vector  fields 

Suppose  that  experimental  3D  image  data  of  compressed  fiber 
systems  for  compression  rates  ce  (0, 10, 20, 30}  (in  %)  are  available, 
see  Figs.  3  and  4.  In  the  following,  we  distinguish  between  two 
kinds  of  compression,  i.e.,  using  a  flat  stamp  and  a  stamp  possessing 

Table  2 

Coefficient  of  determination  R 2  computed  for  the  fitted  vector  fields  of  the 
compression  model  applied  to  experimental  non-woven  GDL  data  for  different  rates 
of  compression  c  and  stamp  geometries. 

c  =  10  c  =  20  c  =  30 

Flat  stamp  0.93  0.97  0.96 

Flow-field  stamp  0.95  0.96  0.97 
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Compression  under  Compression  under 

land  of  flow-field  channel  of  flow-field 


Fig.  11.  Blue  points:  estimated  transition  vectors  by  the  fitting  technique  introduced  in 
Section  3.2  and  applied  to  the  10%  compressed  GDL  by  a  flow-field  stamp.  Surface: 
linear  interpolation  between  these  points.  (For  interpretation  of  the  references  to  color 
in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 

a  flow-field  structure,  see  Fig.  1.  For  each  type  of  stamp,  we  sepa¬ 
rately  apply  the  fitting  technique  introduced  in  Section  3.4  in  order 
to  describe  the  microstructure  of  GDL  compressed  by  an  arbitrary 
rate  of  compression.  Recall  that  the  uncompressed,  experimental 
GDL  is  given  in  a  domain  size  1314  x  1314  x  200  pm3  for  the  series 
of  images  which  are  compressed  by  a  flat  stamp,  and 
1314  x  1314  x  210  pm3  for  the  flow-field  stamp. 

4.1.1.  Flat  stamp 

In  case  of  the  flat  stamp,  we  choose  n  =  1  in  (Section  3.4),  and 

(c))(x,y,z)  =  Pi(c)  z.  Table  2  illustrates  that  the  choice  of  this 
parametric  form  is  reasonable  since  the  coefficients  of  determina¬ 
tion  R2  are  very  close  to  1.  Moreover,  it  turns  out  that  the  function 
P/(c)  =  1  -  c  is  a  suitable  choice,  see  Fig.  12  (left).  Thus,  for  the  flat 
stamp,  {Vc,  ce  [0,  cmax)}  is  given  by  Vc(x,y,z)  =  (1  -  c)z  for 
ce  [0,  cmax).  In  other  words,  by  compressing  GDL  using  a  flat  stamp, 
fibers  are  translated  linearly  according  to  their  locations  with 
respect  to  the  z-coordinate.  Note  that  in  fact,  this  behavior  is 
intuitive  but  it  was  never  shown  before  and  often  used  as  model 
assumption,  see  e.g.  Ref.  [20]. 

4.1.2.  Flow -field  stamp 

In  the  next  step,  we  consider  the  estimation  of  the  parametrized 
vector  field  Vc  for  compression  by  a  flow-field  stamp.  First  note  that 


the  channel  of  the  flow-field  stamp  has  a  width  of  w  =  790  pm. 
Since  the  profile  of  the  flow-field  stamp  replicates  in  y-direction, 
we  assume  that  the  vector  field  Vc  for  compression  by  a  flow-field 
stamp  does  not  depend  on  the  y-coordinate,  but  on  the  orthog¬ 
onal  distance  d(-)  (given  in  pm)  to  the  center  of  the  flow-field 
channel  (orientated  in  (0, 1,0) T  direction  as  displayed  in  Fig.  4). 
Note  that  the  estimated  vector  field  contains  all  information  to 
describe  the  transition  of  fibers  caused  by  the  compression.  Hence, 
all  structural  changes  of  microstructures  that  occur  by  the  uniaxial 
compression  can  be  explained  by  means  of  the  vector  field.  To  get 
an  idea  for  the  choice  of  the  parametrized  vector  field  Vc  concerning 
compression  by  a  flow-field  stamp,  consider  Fig.  11.  In  this  figure, 
we  see  that  fibers  underneath  the  rips  of  the  flow-field  (i.e.  fibers 
with  d(x)  >  w/2)  are  translated  linearly  according  to  their  locations 
with  respect  to  the  z-coordinate.  But  fibers  underneath  the  channel 
of  the  flow-field  (i.e.  fibers  with  d(x)  <  w/2)  follow  a  non-linear 
translation  rule.  Note  that  values  of  Vc  smaller  than  zero  indicate 
decompression  of  the  fibers.  Thus,  fibers  located  close  to  the  center 
of  the  flow- field  channel  are  decompressed,  see  Fig.  11.  Moreover, 
we  observe  that  the  decompression  of  fibers  within  the  flow-field 
channel  is  smaller  (larger)  the  closer  the  fibers  are  located  to  the 
bottom  (top)  of  the  GDL.  Thus,  in  this  case  the  continuous  vector 
field  Vc  is  constructed  by  distinguishing  between  the  two  disjoint 
regions  and  describe  the  translation  of  fibers  in  those  regions, 
separately.  To  assure  the  continuity  of  Vc,  we  include  a  small  tran¬ 
sition  region  (d(x)e  [w/2  -  4,  w/2])  connecting  regions  1  and  2, 
linearly.  Thus,  for  the  parametric  vector  field  Vc,  it  seems  reasonable 
(see  Table  2)  to  put  n  =  6  and  consider  the  set-up 

Vf'(c)'-'Mc)\x,y,z)  =g1(z;c)lld(x)>»+  (gi(z;c) 

+p±^zg^).(<w-S) 

1 1^-20  <d(x)<%  +&2  C)  1  Id(x)<f-20> 

where 

g\(z-c)=pA  (c)-z, 

g2  (d(x),z;c) =p2  ( C)z+p3  (c)d(x)z+p4(c)z2  +p  5  (c)z3  +p6(c)d(x)2z2 

and  ce  [0,cmax).  Moreover,  we  putpi(c)  =  1  -  c  and  for  ie{2,  ...,6} 
we  use  the  family  of  equations  pfc)  =  a1(exp(bIc)  -  exp(d;c)),  where 
the  estimated  coefficients  are  specified  in  Table  1.  Table  2  and  Fig.  12 
(right)  illustrate  that  the  choice  of  Vc  seems  to  be  adequate. 


rate  of  compression  c  in  % 


rate  of  compression  c  in  % 


Fig.  12.  Fitted  functions  p,(c)  for  the  compression  by  a  flat  stamp  (left)  and  a  flow-field  stamp  (right),  where  p2,  ....  P6  are  appropriate  scaled  in  order  to  check  the  goodness  of  fit. 
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Fig.  13.  First  row:  experimental  3D  image  data  of  non-woven  GDL  compressed  by  flat  stamp  (left)  and  by  a  flow-field  stamp  (right)  with  compression  rate  c  =  30%;  second  row: 
extracted  fiber  system  virtually  compressed  by  a  flat  stamp  (left)  and  by  a  flow-field  stamp  (right)  using  the  compression  model  with  compression  rate  c  =  30%. 


4.2.  Validation  of  compression  model 

In  Section  4.1,  vector  fields  of  the  compression  model  for 
compression  by  a  flat  stamp  and  a  stamp  with  flow-field  structure 
are  estimated  in  order  to  describe  compressed  microstructures  of 
non-woven  GDL  in  PEMFC  for  arbitrary  compression  rates 
ce  [0,  cmax).  We  now  check  if  the  compression  model  reproduces  the 
microstructures  of  compressed  non-woven  GDL  for  the  two  stamp 
types  sufficiently  well.  More  precisely,  to  generate  the  virtually 
compressed  microstructures  we  first  consider  the  extracted  fiber 
system  of  the  uncompressed  GDL,  see  Section  2.3.  Subsequently, 
the  extracted  fibers  are  translated  by  the  compression  model  with 
the  help  of  the  estimated  vector  fields  for  the  compression  rates 
ce  {0, 10, 20, 30}  (in  %),  computed  in  Section  4.1. 

From  Fig.  13  it  can  be  seen  that  the  virtually  compressed  images 
coincide  quite  nicely  with  the  synchrotron  images  of  experimen¬ 
tally  compressed  materials.  In  addition  to  visual  comparison  we 
check  the  goodness-of-fit  more  formally.  That  means,  the 
compression  model  is  validated  by  comparing  structural  charac¬ 
teristics  computed  for  experimentally  compressed  and  virtually 
compressed  microstructures. 

In  particular,  for  each  compression  rate  c,  we  compare  the 
porosity  and  the  mean  spherical  contact  distance  for  experimentally 
compressed  and  virtually  compressed  microstructures.  Note  that 
the  porosity  of  a  3D  microstructure  given  in  an  observation  window 
W  is  defined  as  the  volume  of  the  pore  phase  divided  by  the  volume 
of  W  Moreover,  the  mean  spherical  contact  distance  is  the  mean 
value  of  the  closest  distances  from  all  points  located  in  the  pore 
phase  to  the  fiber  phase.  When  describing  the  properties  of  porous 
materials  the  spherical  contact  distance  is  of  high  significance. 


Fig.  14.  The  blue  boxes  display  the  decomposition  in  two  quite  homogeneous  sub- 
regions,  i.e.,  a  channel  and  a  land  region.  (For  interpretation  of  the  references  to  co¬ 
lor  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


Since  it  describes  the  spatial  (micro-)  structure  of  the  pore  phase 
this  information  is  closely  related  to  the  pore  size  distribution.  The 
influence  of  compression  on  pore  size  distribution  and  the  possible 
consequences  for  mass  transport  are  discussed  in  more  detail  in  the 
recent  paper  [23].  Since  the  microstructures  compressed  by  a  flat 
stamp  are  homogeneous,  the  porosity  and  the  mean  spherical 
contact  distance  are  computed  in  the  whole  observation  window. 
The  microstructures  compressed  by  a  stamp  with  flow-field 
structure  can  be  decomposed  in  two  quite  homogeneous  sub- 
regions:  the  region  underneath  the  flow- field  channel  and  the  re¬ 
gion  underneath  the  rips  of  the  flow-field,  see  the  blue  boxes  in 
Fig.  14.  Thus,  for  each  sub-region  separately,  the  porosity  and  the 
mean  spherical  contact  distance  are  computed.  The  results  of 
comparing  structural  characteristics  for  experimentally  com¬ 
pressed  and  virtually  compressed  microstructures  can  be  seen  in 
Fig.  15.  For  both  types  of  stamps  a  good  agreement  is  found.  Note 
that  this  also  justifies  the  model  assumption  of  the  compression 
model,  i.e.,  that  fibers  are  only  translated  in  z-direction  (the  di¬ 
rection  of  the  uniaxial  compression).  Re-considering  Fig.  15  we 
clearly  see  that  the  porosity  (mean  spherical  contact  distance) 
computed  for  the  samples  compressed  by  a  flat  stamp  and  under¬ 
neath  the  land  regions  when  referring  to  compression  by  a  flow- 
field  stamp  decreases  (increases)  linearly  with  respect  to  the 
compression  rate.  In  contrast,  the  porosity  (mean  spherical  contact 
distance)  underneath  the  channel  regions  of  samples  compressed 
by  a  flow-field  stamp  shows  a  non-linear  behavior.  This  can  be 
explained  by  the  computed  vector  fields  considered  in  the  previous 
section  where  underneath  the  land  (channel)  of  the  flow-field  a 
linear  (non-linear)  translation  of  fibers  is  observed.  Note  that  the 
focus  of  the  present  paper  is  on  the  introduction  of  a  compression 
model  which  is  capable  to  generate  microstructures  of  GDL  with  an 
arbitrary  degree  of  compression  in  very  large  volumes  taking  into 
account  the  influence  of  the  flow-field  geometry.  Changes  in  the 
microstructure  that  result  from  compression  and  possible  conse¬ 
quences  for  the  fuel  cell  operation  are  discussed  in  more  detail  in 
the  recent  paper  23]. 

5.  Application  to  simulated  data 

In  addition  to  the  scenarios  considered  in  Section  4,  the 
compression  model  can  be  applied  not  only  to  experimental  data 
but  also  to  virtual  non-woven  GDL  data  sampled  from  stochastic 
microstructure  models.  The  main  advantages  of  virtual  GDL  are  that 
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Fig.  15.  Validation  of  compression  model  by  comparing  the  porosities  (top)  and  the  mean  spherical  contact  distances  (bottom)  computed  for  the  experimentally  compressed  image 
data  (black  lines),  for  data  generated  by  virtually  compressing  the  extracted  fiber  systems  (red  lines)  and  for  data  gained  by  compressing  virtual  non-woven  GDL  by  the  compression 
model  (blue  lines),  where  we  distinguish  between  compression  by  a  flat  stamp  (left)  and  a  flow-field  stamp  (right).  (For  interpretation  of  the  references  to  color  in  this  figure  legend, 
the  reader  is  referred  to  the  web  version  of  this  article.) 


Fig.  16.  First  row:  experimental  (left)  and  synthetic  (right)  uncompressed  non-woven  GDL  data;  second  row:  non-woven  GDL  data  compressed  by  a  flat  stamp  (with  compression 
rate  c  =  30%)  experimentally  (left)  and  using  the  compression  model  (right);  third  row:  non-woven  GDL  data  compressed  by  a  flow-field  stamp  (with  compression  rate  c  =  30%) 
experimentally  (left)  and  applying  the  compression  model  on  a  virtual  fiber  system  sampled  from  the  stochastic  non-woven  GDL  model  (right). 
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they  can  be  produced  in  arbitrarily  large  observation  windows  with 
arbitrarily  many  replications  and  that  their  generation  is  cost-  and 
time-saving.  Moreover,  by  systematically  varying  the  parameters  of 
the  stochastic  microstructure  model,  new  virtual  GDL  can  be  pro¬ 
duced  and  in  combination  with  computational  transport  simula¬ 
tions,  microstructure  and  functionality  relationships  of  the  GDL  can 
be  investigated.  Thus,  after  the  estimation  of  the  vector  fields  as 
discussed  in  Section  4.1,  the  compression  models  for  flat  stamps 
and  stamps  with  flow-field  structure  are  applied  to  (uncom¬ 
pressed)  virtual  non-woven  GDL.  For  the  generation  of  the  (un¬ 
compressed)  virtual  non-woven  GDL  a  stochastic  simulation  model 
is  used  which  has  recently  been  developed  in  Ref.  [22].  In  particular, 
it  has  been  shown  in  Ref.  [22]  that  transport-relevant  characteris¬ 
tics  are  captured  by  the  microstructure  model  quite  well.  Thus, 
simulated  3D  microstructures  drawn  from  the  model  proposed  in 
Ref.  [22  can  be  seen  as  realistic  representations  of  experimental 
uncompressed  GDL.  The  combination  of  this  stochastic  simulation 
model  with  the  compression  model  proposed  in  the  present  paper 
can  be  interpreted  as  a  stochastic  4D  model  which  describes 
continuously  compressed  3D  microstructures  of  non-woven  GDL  in 
PEMFC  for  arbitrary  compression  rates  ce  [0,  cmax). 

By  comparing  the  same  structural  characteristics  as  considered 
in  Section  4  we  validate  the  4D  model,  i.e.,  we  check  if  the  sto¬ 
chastic  4D  compression  model  describes  experimentally  com¬ 
pressed  microstructures  of  non-woven  GDL  adequately.  As  shown 
in  Fig.  15,  a  very  good  accordance  between  these  characteristics  is 
found.  Additionally,  there  is  a  nice  visual  accordance  of  experi¬ 
mental  3D  image  data  of  compressed  non-woven  GDL  and  the 
virtually  compressed  non-woven  GDL  sampled  from  the  stochastic 
simulation  model  proposed  in  Ref.  [22]  in  combination  with  the 
compression  model  introduced  in  the  present  paper,  see  Fig.  16. 

6.  Conclusions 

We  have  introduced  a  novel  parametrized  model  that  describes 
the  3D  microstructure  of  compressed  fiber-based  materials.  Given 
the  input  of  a  3D  microstructure  of  some  fiber-system,  the  model 
compresses  the  system  of  fibers  in  a  uniaxial  direction  for  arbitrary 
compression  rates.  The  compression  model  was  applied  to  non- 
woven  GDL  in  PEMFC  and  subsequently  validated  by  comparing 
structural  characteristics  computed  for  experimentally  compressed 
and  virtually  compressed  microstructures,  where  two  kinds  of 
compression  -  using  a  flat  stamp  and  a  stamp  with  flow-field 
profile  -  were  considered.  For  both  types  of  stamps,  an  excellent 
agreement  between  experimental  compression  and  virtual 
compression  was  found.  In  addition,  the  compression  model  was 
combined  with  a  stochastic  model  for  uncompressed  non-woven 
GDL.  This  allows  to  generate  virtual  GDL  in  arbitrarily  large  vol¬ 
umes  with  arbitrarily  many  replications  and  with  low  computa¬ 
tional  efforts.  Thus,  by  systematically  varying  the  parameters  of 
these  combined  models,  new  virtual  GDL  can  be  generated  and  in 
combination  with  computational  transport  simulations,  the 
microstructure-functionality  relationship  of  the  GDL  can  be  inves¬ 
tigated.  Currently,  the  compression  model  is  used  for  virtual  sce¬ 
nario  analyzes  with  the  general  aim  to  quantify  the  effect  of 
compression  using  a  stamp  with  a  flow-field  structure  compared  to 
compression  by  a  flat  stamp.  The  results  of  this  ongoing  research 
will  be  reported  in  a  forthcoming  paper.  But  note  that  the  presented 


compression  model  is  very  flexible  and  can  thus  be  applied  to 
various  problems  not  only  in  the  field  of  GDL  in  PEMFC. 
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